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ABSTRACT 

One-armed oscillation modes in the circumstellar discs of Be stars may explain the 
cyclical variations in their emission lines. We show that a three-dimensional effect, 
involving vertical motion and neglected in previous treatments, profoundly influences 
the dynamics. Using a secular theory of eccentric discs that reduces the problem to a 
second-order differential equation, we show that confined prograde modes are obtained 
for all reasonable disc temperatures and stellar rotation rates. We confirm these results 
using a numerical analysis of the full set of linearized equations for three-dimensional 
isothermal discs including viscous terms that couple the horizontal motions at different 
altitudes. In order to make these modes grow, viscous damping must be overcome by 
an excitation mechanism such as viscous overstability. 



Key words: accretion, accretion discs 
stars: emission-line. Be 



circumstellar matter — hydrodynamics - 



1 INTRODUCTION 



Classical Be stars (jPorter fc Rivinius|[2003l ) are rapidly ro- 
tating early-type stars that exhibit Balmer emission lines. 
It is widely agreed that these lines, which are generally 
double-peaked, originate in a relatively thin circ umstellar 
disc that is in approximately Keplerian rotation (|Okazakil 
I2OO7I . and references therein). While the precise mechanism 
by which the disc is formed remains controversial, it is likely 
to resemble a viscous decretion disc that is expelled by the 
action of a torque at i ts inner boundary l|Lee et al.lll99ll : 
IPorter fc Riviniusll2003l ). 

Many Be stars show cyclical variations in their double- 
peaked emission lines over years or decades, with th e red and 
blue peaks alternately becoming more prominent jOkazakil 
I1997I . and references t herein). An exp lanation of this phe- 
nomenon was given by lOkazakil (|l99ll ), who p roposed that 
a low-frequency, one-armed oscillation mode l|Katol Il983l ) 
occurs in the disc. This is equivalently to saying that the 
disc becomes eccentric and the slow precession of its ellipti- 
cal shape gives rise to the cyclical changes in the observed 
emission lines. 

Eccentric discs, in which fluid elements, solid particles 
or stars follow elliptical orbits of variable eccentricity around 
a central mass, have further applications in systems as di- 
verse as planetary rings, protoplanetary systems, close bi- 
nary stars and galactic nuclei. A detailed understanding of 
the origin of eccentricity and the rates of precession in the 
circumstellar discs of Be stars would therefore be of general 
interest. 

lOkazakil originally considered a disc that orbits 



in a point-mass potential and obtained a sequence of ret- 
rograde modes in which the precession of the disc is in a 
direction opposite to its rotation. Retrograde precession is 
a natural consequence of the pressure forces in the disc, 
which cause a small departure from Keplerian rotation and 
allow the eccentricity to p ropagate in a wavelike manner. 
The global modes found bv lOkazakil l|l99lh are weighted to- 
wards the outer part of the disc and their periods become 
extremely long as the outer radius of the disc is increased to 
realistic values. 



[ Papaloizou et al.l (|l992|) and ISavoniie fc HeemskerkI 
(|i99a r considered the effect of the quadrupole gravitational 
potential associated with the rotational deformation of the 
star, which tends to cause a prograde precession of ellipti- 
cal orbits. They showed that, when the quadrupole effect is 
taken into account, prograde modes can be obtained that 
are naturally confined in the inner part of the disc and are 
insensitive to the outer boundary condition. Subsequent ob- 
servations confirmed that the precession is indeed prograde 
(TcltinE ct a], , 1994). 

However, lOkazakil 11993) found that confined prograde 
modes can be obtained only when the disc is sufficiently cool, 
so that the quadrupole effect dominates over the tendency 
of pressure to produce retrograde precession and extended 
modes. He concluded that a different mechanism is required 
in the hotter discs of early-type Be stars such as the pro- 
totype 7 Gas, and proposed that radiative line forces could 
explain the required prograde precession. Unfortunately the 
modelling of radiative force s is su bject to considerable un- 
certainty. iFift fc Harmane3 l|2006h found that the resulting 
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model has little predictive power owing to the sensitivity of 

the results to the p arameters. 

More recently, IPaoaloizou fc Savoniid ()2006l ) investi- 
gated an alternative way to obtain prograde modes in hotter 
discs. This involves replacing the rigid inner boundary con- 
dition at the stellar surface with a free boundary condition, 
on the basis that a gap is formed between the star and the 
disc. Such a gap is not generally expected in the scenario of 
the viscous decretion disc but might be possible in alterna- 
tive models. 

All of the treatments described so far are based on two- 
dimensional models of the disc that neglect aspects of its 
vertical structure and motion. At first sight, such an ap- 
proach seems reasonable for studying eccentric modes in thin 
discs, where the motion might be assumed to be purely hor- 
izontal and independent of height. However, in presenting 
a three-dimensional, non-linear theory of eccentric discs, we 
have previously argued th at three-dime nsional effects are of 
considerable importance (jOgilviel I2OOII ). In particular, the 
variation of the vertical gravitational acceleration around 
an elliptical orbit excites an oscillatory vertical motion in 
an eccentric disc that should not be neglected. 

In this paper we show that, in fact, three-dimensional 
effects are essential to understanding the precession of ec- 
centric discs around Be stars. They allow confined prograde 
modes to be obtained even when the stellar quadrupole mo- 
ment is negligible and when the inner boundary is rigid. This 
property allows us to give a unified description of eccen- 
tric discs of Be stars of all stellar types without introducing 
uncertain radiative forces or modifying the inner boundary 
condition. While we do not deny that in some cases radia- 
tive forces might be important, or that the inner boundary 
condition might differ from a rigid one, we show that these 
innovations are unnecessary. 

Most previous analyses have not discussed the pro- 
cesses that could cause eccentric modes to grow or decay, 
and which are therefore relevant to explaining the occur- 
ren ce of eccentric discs around Be stars . Viscous overstabil- 
ity l|Katdll97 8': 'Latter fc Ogilvic"2006l') provides a possible 
explanation; [Nogucrucla et al. (2001) have applied this idea 
to Be stars and estimated the associated growth rate. We 
defer to a future investigation a detailed analysis of the ef- 
fects of viscous or turbulent stresses. A non-linear treatment, 
which could address the saturation of the growth mecha- 
nism and attempt to predict the observed amplitudes and 
precession rates of the eccentric modes, also remains to be 
undertaken. 

The structure of this paper is as follows. In Section [2] 
we describe the basic state of the disc. We then discuss an 
approximate, secular theory of three-dimensional eccentric 
discs in Section [3l comparing it with the two-dimensional 
theories used by other authors. In Section |4] we solve the 
full linearized equations accurately using a spectral method, 
including some effects of viscosity, and compare the results 
with those of the secular theory. Conclusions are given in 
Section [31 



2 BASIC STATE OF THE DISC 

We consider a basic state consisting of a steady, axisymmet- 
ric disc around a rotating star. Adopting cylindrical polar 



coordinates {r,(f),z), we write the gravitational potential of 
the star as 
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where M and R are its mass and (equatorial) radius. This 
potential consists of monopole and quadrupole components, 
the latter being of dimensionless strength Q and arising from 
the rotational deformation of the star. (The parameter Q is 
related to the gravitational moment J2 used in planetary 
science by Q = 3J2/2.) For a star with uniform angular 
velocity f2,, we have Q = k2^'iR^ /GM, where is the 
apsidal motion constant. We neglect higher-order multipole 
components as well as the self-gravitation of the disc. 

In common with most other treatments, we assume that 
the stellar radiation maintains the disc at a constant tem- 
perature Td, slightly less than the effective temperature Tc 
of the star. The isothermal sound speed Cs — {TZTd/ /-t)^^^ is 
therefore constant, both in the equilibrium state and for the 
perturbations. We neglect any viscous or turbulent stresses 
and any resulting meridional motions in the disc. 

The basic state then has velocity u = rf2(r-)e0. The 
angular velocity Q is independent of z because the basic 
state is isothermal and therefore barotropic. The balance of 
forces requires rQ^ = V($ -f /i), where /i = Cg Inp is the 
(pseudo-)enthalpy of an isothermal gas. 

The potential and enthalpy in the midplane are 
-I>m(r) = -l>(r,0) = -{GM/r)[l + {Q/3){R/rf] and hn.{r) = 
h{r,0), and the radial force balance implies 



GM 
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The associated epicyclic frequency K(r) is given by 
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We also refer to the Keplerian angular velocity QK{r) defined 
by 

«i - ^. (4) 

For a thin disc we may expand the potential about the 
midplane to obtain $ ~ $m + ^fi^z^, where the vertical 
frequency Qz{r) is given by 



GM 



1 + 



3«(f) 



(5) 



The vertical force balance then implies h = hm — ^fi^z^ and 
so p = exp[~{z^ /2H^)], where Pmir) is the midplane 
density and H{r) — Cs/fiz is the scaleheight. The surface 
density is E(r) = {2-kY^'^p^H. 

Since the basic state is steady and axisymmetric, wave 
modes may be considered in which the dependence on az- 
imuth and time is of the form eiqy{im(f) — iij-it), where m is 
the azimuthal wavenumber and u) is the wave frequency. In 
the case m = 1 that is of interest here, uj is also the angular 
pattern speed of the mode in an inertial frame of reference. 
This can be identified with the precession rate of the eccen- 
tric disc, which is positive if the precession is prograde (i.e. 
in the same direction as the rotation of the disc). 
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3 SECULAR TREATMENT 

3.1 Comparison of two- and three-dimensional 
theories 

An approximate theory can be derived in which the disc is 
assumed to be nearly Keplerian and the frequency uj of the 
wave is much less than the orbital frequency il of the disc. 
This type of approximation, which is related to the secu- 
lar theory of celestial mechani cs, has be e n discussed in the 
two-dimensional linear case biJ^emai^ (|200l),|PaEaioiio3 
l|2002t) and iGoodchild fc Ogilvid hOO&i. and in the three- 
dimensional non-linear case VlPEilvic (2003). The deriva- 
tion of this theory in the case of a three-dimensional isother- 
mal disc around a Be star is presented in Appendix [X] 
which draws on the analysis of Section |3] below. The two- 
dimensional approximation is also discussed there. 

In the secular theory the radial velocity u{r) satisfies 
the equation 



{UJ - f)u = —dr 



where /(r) is given by 



-d, 



20 

in a two-dimensional disc, but 
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in a three-dimensional disc. With the dependence e 



(6) 
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assumed above, the radial velocity is related to the com- 
plex eccentricity E{r) — ee'", where e(r) is the eccentricity 
and m{r) is the longitude of periastron measured in a frame 
of reference that rotates with the angular pattern speed ui 
of the mode, by u* = ir^lE (Ogilvic 2001). It can be seen 
from equation © that / is a local contribution to the global 
precession rate u> of the eccentric mode. Indeed, an integral 
expression for the mode frequency in the case of rigid bound- 
ary condition^ (m = at r = rin and r = rout) is 



T,r\u[ 



dr = 
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dr 



(9) 



which follows from equation (|6]) after multiplication by 
Ert6*/0 and an integration by parts. The first term on the 
right-hand side shows the contribution to the precession rate 
from / in the form of an integral weighted by the structure 
of the mode, while the second term shows a retrograde con- 
tribution associated with pressure. (Note that the pressure 
also contributes indirectly through its effect on /.) 

In a two-dimensional disc / is given by equation (O and 
corresponds to the expected expression for the local apsidal 
precession rate. Note that, in this case, / « — k if , as as- 
sumed, the disc is nearly Keplerian — <C ^l^). There 



^ Altliougli there is no physical justification for a rigid outer 
boundary condition, a confined mode decays suflBciently fast with 
increasing r that the same integral relation is recovered in the 
limit of large rout. Such modes are in fact insensitive to the outer 
boundary condition. 



is, in f act, more than one ve r sion o f the two-dimensional 
theorv. lSavoniie fc HeemskerkI l|l993h work throughout with 
vertically averaged equations and find 
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-9r(r dr In E). 
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lQkazakill|l99'i1 ) calculates D, and k using a three-dimensional 
equilibrium disc but then applies vertically averaged equa- 
tions for the perturbations. He therefore uses 



f = Q 



GMW 



-dr{r 9rlnpm), 



(11) 



which agrees with equations ([2} and ((3} above. 

In a three-dimensional disc, however, equation ([8} gives 
the relevant expression for / as 



f = Q 



GMW 



c 9c 
-dr{r^dr Inpni) + ■ 



(12) 



The last term represents an additional local contribution 
to the prograde precession of the disc that arises from the 
three-dimensional dynamics including the vertical motion; 
it is discussed further in Section [3.71 below. For the param- 
eters relevant to Be stars, this term is never negligible (it 
corresponds to a precession period of order 1 yr at the stel- 
lar surface) and declines much more slowly with radius than 
the quadrupole term. 

We refer to the three theories described above as 
2DS' (equation \T0 \ above; ISavoniie fc Heemsker^ Il993l : 



Paoaloizou fc Savoniid l2006l i . '2DO' (equation [TT] above 



Okazakilll99ll . ll997h and '3D' (equation [l2l. Note, however. 



that the secular approximations were usually not employed 
in the cited papers. The 3D secular theory is based on similar 
as sumptions and approximations to the non-linear analysis 
of lOgilvid l|200ll ). In particular, some viscous or other stress 
is required to couple different layers in the disc so that they 
tend to adopt the same eccentricity. In the absence of such 
stresses the eccentricity may vary significantly with z and 
the results can differ. This complication is discussed in Sec- 
tion |4l where the full linearized equations are solved. 



3.2 Application to power-law discs 

It is consistent with the spirit of the secular approximation 
to neglect the differences between il, k, Q^, and JIk except 



where essential (i.e. in the quantity O 



Accordingly, 



we may take D, oc r~^^^ and H oc r^^^. For a surface density 
profile E oc r""^ the density in the midplane varies as pm oc 
J..-CT-3/2 -yyr^ introduce the dimensionless radial coordinate 

^ = J (13) 
and the small parameter 
/ R 

which is a measure of the angular semithickness H/r of the 
disc at the stellar surface. The three theories then give 

/ 

n 
I 



(14) 
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(18) 



We therefore consider the generic form 
/ ^ 

n 

where s is a constant. It is convenient to work with di- 
mensionless, rescaled values of the mode frequency and 
quadrupole strength, uj and Q, defined by 



Equation ((Gjl then becomes 



(19) 
(20) 



-7/2 



+ SX 



[x-+'dAx--+''\)\ 



(21) 



[Cj - {Qx- 
1 

which is an eigenvalue problem of Sturm-Liouville form, 
when considered together with appropriate boundary condi- 
tions. Note that the parameter e drops out of the equation 
under the above rescalings. For prograde modes (iIj > 0), one 
solution of this equation decays exponentially as a; oo, 
while the other solution grows exponentially. We therefore 
consider the equation in the domain 1 < a; < oo, requiring 
solutions to satisfy a rigid boundary condition u = at the 
stellar surface x = 1 and to decay as a; ^ oo. 

Reasonable v alues of Q for Be stars can be estimated 
foUowing Fir t fc Harmaned (120061 ). Their Table 4 is based on 
stellar models of intermediate main-sequence age. There is 
considerable uncertainty in the applicability of these models 
because of the rapid rotation of Be stars. Nevertheless, using 
these data, we estimate that Q ranges from approximately 8 
at the lower end [M = 2.51 Mq) to approximately 12 at the 
upper end (A/ = 15.85 Mq). These values assume, somewhat 
arbitrarily, a disc temperature of Td = (2/3)re and a stellar 
rotation rate that is 95% of the critical value. Similarly, for 
the stell ar models of 7 Cas and 5 9 Cyg described in Tables 2 
and 3 of lFift fc Harmaned (|2006l ). we estimate Q ~ 9 and 13, 
respectively. It is clear that significantly larger values of Q 
cannot be obtained by further increasing the stellar rotation 
rate, and in fact smaller values may be more appropriate. 
The values of e for all these models are in the range 0.021- 
0.026, again assuming that Td = (2/3)rc. 



3.3 Solutions in the absence of a quadrupole 

Consider first the case without a quadrupole term, Q = 0. 
The solution that decays as a:: 00 is then 

u = x^''-''^'^K^{z), 



with 



u = 2[(a + 2f - Ssf'\ z = A{2Gjf/^x^'\ 



(22) 



(23) 



where i s the modified Bessel functio n of the second kind 
of order v l|Abramowitz fc Steguniri965l ). Note that z is real 
and positive for the prograde modes of interest. This solution 
has no zeros in z > if is real, but does so if u is imaginary. 
To match a rigid boundary condition at a; = 1 we therefore 
require < 0. The 2DS theory gives 8s — Aa, so = 
4((T^ + 4) and modes are never confined. The 2DO theory 
gives 8s = 4(7 + 6, so = 4(cr^ — 2) and modes may be 
confined for cr^ < 2. The 3D theory gives 8s = 4(t -I- 24, so 
= 4((T^ — 20) and modes may be confined for a'^ < 20. 
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Figure 1. Confinement of prograde modes in the absence of a 
quadrupole. The modified Bessel function K,^{z) (multiplied by 
a constant for convenient normalization) is plotted against x = 
r/R for the cases v = 81 (solid line) and = VSi (dashed line). 
According to equation II22I I. the modified Bessel function, which 
decays exponentially as a: —> 00 if f is imaginary, describes the 
confinement of the eigenfunction. The solid line is relevant to 
the 3D theory for a surface density index tr = 2, and shows a 
good confinement of the mode. The dashed line is relevant to the 
2DO theory in the (optimal) case cr = 0, but here the confinement 
is very poor. 



This analysis is somewhat misleading because the 
modes may not be adequately confined to be applicable to 
the discs of Be stars. In Fig. [1] we plot K,^ (z) versus x for 
= 8i and 1/ — \/8i, choosing Cj such that K,y{z) = at 
X — 1. Only the former case provides an adequately con- 
fined mode. This shows that the 2DO theory cannot in 
fact produce confined prograde modes in the absence of a 
quadrupole, even for the most optimistic choice of surface 
density profile, because {ul is never large enough. However, 
sufficiently large imaginary values of 1/ are obtained in the 
3D theory. 

When f is imaginary, confined solutions are obtained, 
in principle, when Cj = 2^/32 where Zn.i^ is the nth zero 
of K,^{z) in 2 > 0. For example, consider the 3D theory 
with a — 2, which is the value expected for a steady decre- 
tion disc with constant alpha viscosity parameter, far inside 
its outer radius. Then s = 4 and v — 8i, and zeros oc- 
cur a.t z = 4.802, 3.067, 2.029, . . . Therefore, in principle, 
Q = 0.7207, 0.2940, 0.1286, ... are the scaled frequencies 
of confined modes. However, only the first of these, corre- 
sponding to an eigenfunction with no nodes (Fig. [1] solid 
line), gives rise to a mode that is adequately confined in a 
disc of modest radial extent. 



3.4 Critical quadrupole strength 

The integral expression for the dimensionless frequency 
eigenvalue of a mode that satisfies a rigid boundary con- 
dition at the stellar surface and decays as r 00 is (cf. 



Three-dimensional eccentric discs around Be stars 5 



equation [9)l 

^-CT+5/2^2 _ / (Qx"'^""'" + SX~'^^^)u^ dx 

x^+'^id^x-'^'/^)]' dx, (24) 

where we take u to be real. FoUowing lPapaloizou fc Savoniiel 
l|2006t ). we note that this expression has the usual variational 
property associated with self-adjoint eigenvalue problems. In 
this confined prograde mode exists if and only the 

right-hand side of this equation can be made positive by 
a trial function u{x) satisfying the appropriate boundary 
conditions. This is clearly possible if either Q or s is large 
enough. What is the minimum value of Q (for a given value 
of s) such that the right-hand side can just be made to vanish 
for a non-trivial u7 This happens when 

x~'^~^u^dx = - J x'^'^^ ^dx{x~'^^^'"^u)'\ dx 

/oo 
sx-'^+'^u^ dx. (25) 

The minimum value of Q for which this equation can be sat- 
isfied is given by the corresponding Euler-Lagrange equa- 
tion, which is just equation (|2ip with Cij set to zero. The 
solutions are 

u = x^^-^^^^J±^{w), (26) 
with 

^^6=3'^" + ^) ' "^ — (2; -(2^^ 

If ^2 < we have seen that confined modes can be found, 
in principle, even if Q = 0. Consider then the case /i^ > 0. 
The solution that decayfl as x — > oo is J^{w). The critical 
condition for a rigid inner boundary condition is therefore 
Q = 9wft/8, where is the first zero of J^iw) in w > 0. 
(This is an increasing function of /i and therefore a decreas- 
ing function of s.) 

For example, if ct — 2, we require Q > 10.8 for 
confinement in the 2DO theory, or Q > 15.9 in the 
2DS theory. These values are not very sensitive to a; 
IPapaloizou fc Savoniiel l|2006l ) quote Q > 17.3 for a = 5/2 
and a rigid inner boundary condition, which agrees with this 
analysis. 

3.5 Schrodinger analogy 

A useful description of confined modes uses an analogy with 
bound states in quantum mechanics, li y = x^^^ and it — 
^o-/2-i3/8^^j^^^ then '(/'(y) satisfies the Schrodinger equation 

- + {y{y) -E]^p = (28) 
with an effective potential 

2 Since J^(io) oc for small to, m oc x'°'~'^(-'^+^')l''2 for large x. 
This solution satisfies the condition that x''^^udx{x~'^^^'''^u) — ► 
as a; — > oo, which is required to carry out the integration by parts 
and derive the variational principle and Euler-Lagrange equation, 
while J—u does not. 




Figure 2. Effective potentials V{y) (equation I29I I describing the 
confinement of prograde modes by analogy with bound states in 
quantum mechanics. Potentials are shown for the three theories 
for the case <t = 2 and Q = 10. The bottom of the potential well 
at y = 1 (i.e. at the stellar surface) in each case is considerably 
deeper than can be shown on the graph. However, only in the 
3D theory is the well deep and wide enough in this case to support 
a bound state in the case of a rigid inner boundary condition. 

V = ^(16cr^ + 64(7 + 63 - 128s)y"^ - 32Qy''^^ (29) 

and an effective energy eigenvalue E = — 32d>. The coef- 
ficient of l/4y^ in V is (16(t^ -f- 63) for the 2DS theory, 
(16(7^ - 33) for 2DO and (16cr^ - 321) for 3D. It is therefore 
likely to be negative only in the 3D theory. Note that the 
disc occupies the region y > 1. A bound state of negative 
energy, equivalent to a confined prograde mode, can be ob- 
tained if there is a sufficiently deep and wide potential well. 
The Q term tends to create a deep well close to the stellar 
surface. In the 2D theories it competes with the y~^ term, 
which contributes a repulsive potential. In the 3D theory, 
however, the y~^ term creates a much wider well, so al- 
lowing broader modes with slower precession rates, even if 
Q — 0. These potentials are plotted in Fig. [2] for the case 
(7 = 2 and Q = 10. For these parameters the 3D potential 
is deep and wide enough to support a bound state in the 
case of a rigid inner boundary condition, while the others 
are not. 



3.6 Precession rates and mode shapes 

We now compute the eigenvalues Q of equation (|2ip by a 
shooting method, adopting a rigid inner boundary condition 
at the stellar surface and seeking prograde confined modes 
that decay as a; ^ oo. (In practice this is done by imposing a 
rigid outer boundary condition and verifying that the eigen- 
value is completely insensitive to the value of the outer ra- 
dius provided it is sufficiently large.) The results are shown 
in the left panel of Fig. [3] Here we confirm that confined 
prograde modes exist in the 2D theories only for sufficiently 
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Figure 3. Dimcnsionlcss precession rates of confined prograde modes in the three theories, for the case a = 2. The left and right panels 
are for rigid and free inner boundary conditions, respectively. These results are insensitive to the outer boundary condition provided that 
the disc is reasonably large. In the 2D theories, confined prograde modes are found only for sufficiently large values of Q. The dotted 
line indicates a second confined mode, with one radial node in its eigenfunction, found in the 3D theory; this and higher-order modes 
are less well confined and therefore more sensitive to the outer boundary condition. 



large values of Q, as described in Section|3]4] In contrast, the 
3D theory allows such modes to be obtained for any value 
of Q. Furthermore, the precession rates are much larger in 
the 3D theory. Since the realistic values of Q are probably 
in the vicinity of 10 or smaller, it is clear that the three- 
dimensional effects are of essential importance in the case of 
a rigid inner boundary condition. 

In the right panel of Fig. [3] we show comparable re- 
sults for the free inne r boun dary condition considered by 
IPapaloizou fc Savoniid (2006) , meaning that the Lagrangian 
pressure perturbation is zero at r — R. In the secular theory 
th is is equivalent to drE = 0, or u + 2rdrU — 0. As described 
bv lPapaloizou fc Savoniid (|2006l ) , a free inner boundary con- 
dition allows prograde modes to be found for smaller values 
of Q in the 2DS theory. It is still true, however, that the 
three-dimensional effects have an important effect on the 
results. Since the eigenfunctions obtained with a free inner 
boundary condition are generally peaked at the inner radius, 
there may also be a conflict between the nonlinear develop- 
ment of such an eccentric mode and the existence of a stellar 
surface, unless there is a wide gap between the star and the 
disc. 

To convert these eigenvalues into physical units, we 
again make use of the stellar models in iFii't fc Harmane3 
l|2006t ). We then find that the precession period is 



much as the rotation affects the equatorial radius, but Cp 
is inversely proportional to the assumed disc temperature. 

The shapes of the confined modes in the 3D theory are 
illustrated in Fig. 2] For larger values of Q, the mode is 
increasingly confined in the inner part of the disc. Although 
higher-order modes may exist, one of which is referred to 
in Fig. O we focus here on the fundamental confined mode 
with the simplest radial structure. 

Since Q is close to 1 in the 3D theory for reasonable 
values of Q and with a rigid inner boundary condition, pre- 
cession rates in the vicinity of 1-3 yr are obtained. Observed 
cycle times, which can vary significantly ev en for the sam e 
star, are more usually in the range 5-10 yr (|Qkazakilll997l ). 
This discrepancy might occur because the discs have dif- 
ferent temperature or density profiles from those assumed 
here, either of which would affect the precession period. 
Non-linearity of the eccentric mode may also increase the 
precession period by altering the distribution of eccentricity 
so that it is less peaked in the inner part of the disc. Fur- 
thermore, non-isothermal effects or radiative forces might be 
important. In some cases, including 7 Gas and 59 Cyg, the 
presence of a relative close binary companion may also af- 
fect the dynamics of the eccentric mode, although it would 
seem most likely to contribute to the prograde precession 
and therefore to decrease the precession period. 



P=^yr, (30) 

where Cp ranges from approximately 1.6 at the lower end 
to approximately 2.3 at the upper end (or 2.4 for 7 Gas and 
1.8 for 59 Cyg). This conversion factor is independent of 
assumptions regarding the stellar rotation rate, except inas- 



3.7 Physical interpretation of the 
three-dimensional dynamics 

What is the origin of the additional prograde precession 
that occurs in a three-dimensional disc? As described in Ap- 
pendix|X]and Section|4]below, two different types of motion 
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10 20 30 40 50 
r/R 

Figure 4. Eigenfunctions of confined prograde modes in the 
3D secular theory with a rigid inner boundary condition. The 
eccentricity (proportional to u/rQ, and with arbitrary normal- 
ization) is plotted versus radius for the case c = 2, and for rela- 
tive quadrupole strengths Q = (solid line), 5 (dotted line), 10 
(dashed line) and 20 (dot-dashed line). The first case corresponds 
to the solid line in Fig. [Tl 



are involved in an eccentric disc. One (corresponding to the 
n — mode in later sections) consists of horizontal veloci- 
ties and enthalpy perturbations that are independent of z; 
this describes the eccentric orbital motion of the gas. The 
other (corresponding to n = 2) involves a vertical velocity 
proportional to z and an enthalpy perturbation proportional 
to z'^ — H^; this is a vertical 'breathing' mode of the disc. 
Coupling of these motions occurs because of the variation 
of the vertical gravitational acceleration (or, equivalently, 
the vertical frequency flz, or the scaleheight H) with ra- 
dius. Vertical hydrostatic equilibrium cannot be maintained 
in an eccentric disc because a fluid element in an elliptical 
orbit experiences a vertical gravitational acceleration that 
oscillates with the orbital frequency. The breathing mode is 
excited and the associated enthalpy perturbation affects the 
horizontal dynamics, contributing to the precession of the 
eccentric mode. This contribution is found to be always pos- 
itive. In the more general situation of a non-isothermal disc 
undergoing adiabatic perturbations, the three-dimensional 
contribution to / (in which we compare a 3D-type theory 
with a 2DS-type theory) is 



3(7+1) P 
27 Er^n' 



(31) 



where 7 is the adiabatic index and P is the vertically inte- 
grated pressure (Ogilvie & Goodchild, in preparation), but 
the derivation of this expression is beyond the scope of the 
present paper. 



4 FULL TREATMENT 
4.1 Inviscid dynamics 

The dynamical equations governing an isothermal, inviscid 
disc are 



{dt + u-V)/i = V-u. 



(32) 
(33) 



When these are linearized about the basic state described 
in Section [21 we obtain 



— icDuJ, — 2Slu0 — —drfi , 

. , / / imh' 

212 r 

— itbu'^ = —dzh', 

—itbh' + u'^drh + u'zdzh 
2 



1 , 1^^,* / 

-dr(rur) H h dzUz 

r r 



(34) 

(35) 
(36) 

(37) 



where Cj = uj — mQ, is the Doppler-shifted wave frequency 
and the perturbations, denoted by primed quantities, have 
the dependen ce exp(— ia;^ -I- imcf)). 

Following lOkazaki et al.l l| 19871 ). we decompose the ver- 
tical structure of the mode into the basis of Hermite poly- 
nomials defined by 



He„(C) =e' 



C/2 



-CV2 



(38) 



where ( — z/H is a dimensionless vertical coordinate and 
n — 0,1,2,... These polynomials satisfy the differential 
equation 

He:: (C) - C HeUC) + 'I He„ (C) = 0, 
the recurrence relations 
HeUO =nHe„_i(C), 
CHe„(C) =He„+i(C)+nHen-i(C) 
and the orthogonality relation 

e-^'/" He™(C) He„(C) dC = (271)'/"^! 5™„. 



(39) 

(40) 
(41) 

(42) 



The first three are Heo(C) = 1, Hei(C) = C and He2(C) 
We therefore expand 



u'r{r,z) = u„{r) Hen(C), 

n 

u'^{r,z) = ^u„(r)He„(C), 

n 

u'z{r,z) = y^Wn(r) Hen-i(C), 

n 

h'{r,z) = ^ft„(r)He„(C). 



(43) 
(44) 
(45) 
(46) 



with u„ ~ Vn = h„ ~ for 71 < and w„ — for n < 1. 
Bearing in mind that H depends on r, we have 

a.He„(C) = -(a.ln//)[nHe„(C) + n(n-l)He„_2(C)]. (47) 
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The projected equations are then 

— — 20,V„ = —drhn 

+ (9rln//)[n/i„ + (n+ l)(n + 2)/i„+2], 



K irnhn 

lUVn + T^ltn = 

2S2 r 

nhn 

- ILOWn = — , 

£1 



(48) 
(49) 

(50) 



while equation ((51} is unchanged. These terms act to damp 
the mode, but have most effect on components of large n. 
They have no effect on uo and no, which represent horizontal 
motions independent of z. These viscous terms can also be 
thought of as providing a coupling between different layers 
of the disc and thereby encouraging it to adopt a horizontal 
motion independent of z. We show below that this effect is 
of considerable importance. 



ct rS 



+ 



"r If 

0, 



+ idr\nH){nu„+u„-2) = 0, (51) 

(cf. lTanaka et al.ll2002l : IZhang fc Laill2006l ). 

This approach corresponds to a (Galerkin) spectral 
treatment of the partial differential equations governing the 
linearized dynamics, which is much preferable to a finite- 
difference treatment. In practice this system of equations 
must be truncated by setting u„, etc., to zero for n > N 
for some integer A'^. The 2D0 theory is obtained, in fact, by 
considering a radical truncation. A'' = 0, of the equations. 

4.2 Selected viscous effects 

To include viscosity, a term 

iv-T (52) 
P 

should be added to the right-hand side of the equation of 
motion (1321). where 



T = pu[Vu + iVuy] + p(i/b - |i^)(V-u)l 



(53) 



is the viscous stress tensor. In the context of an isothermal 
disc it is reasonable to assume that the kinematic shear and 
bulk viscosities i/ and Vh depend only on r. We parametrize 
them as 



ly = aCsH, 



Vh = ClbCs-ff. 



(54) 



A full treatment of the effects of viscosity is compli- 
cated, not only because the above expression for the viscous 
force must be evaluated in cylindrical polar coordinates and 
then projected on to the basis of Hermite polynomials, but 
also because the basic state is modified to include a merid- 
ional flow driven by viscous stresses, which should be con- 
sidered in the linearized equations. This problem is therefore 
deferred to a future investigation. 

In the present paper we adopt a simpler approach in 
which only selected viscous effects are included. We consider 
what might be assumed to be the dominant viscous terms, 
i.e. those involving two derivatives with respect to z. Since 



1 n 
-a.[pa.He„(C)] = -7HHe„(C), 



(55) 



the inviscid perturbation equations (|48|l - (|50p are modified 
by the addition of the viscous terms 



V 



nun, 



('^b + |i/) 



(n - l)w„, 



(56) 
(57) 

(58) 



4.3 Numerical solutions 

We solve the system of ordinary differential equations in r 
for modes with m — 1 using a Chebyshev collocation (i.e. 
pseudospectral) method. This approach converts the differ- 
ential equations and boundary conditions into an algebraic 
generalized eigenvalue problem for the frequency uj, which 
we solve using a standard direct method. Specifically, equa- 
tions (|48p - (|5ip . supplemented by the viscous terms (|56|l - 
(|58|) . are solved for n = 0,2, i, . . . , N with u„, etc., set to 
zero for n > N. Rigid boundary conditions u„ = are 
adopted at both inner and outer boundaries, but it is en- 
sured that the modes obtained are completely insensitive to 
the value of the outer radius and therefore to the choice of 
outer boundary condition. 

For comparison with the results in Section [3.21 we con- 
sider a disc with a midplane density profile pm oc r~'^~^^'^ . 
We also include a shear viscosity corresponding to a constant 
a parameter, but no bulk viscosity. 

Sample results are shown in Table [T] The convergence 
of the eigenfrequency with increasing truncation order TV of 
the Hermite polynomial basis is remarkable. The case N — 
corresp onds exactly to the two-dimensional theory consid- 
ered bv lQkazakil (Il99ll ). and therefore agrees well with the 
2D0 secular approximation. Here Q = 20 is large enough 
to support a confined prograde mode. The precession rate 
is much larger in the case N = 2 and hardly varies as fur- 
ther Hermite polynomials are included. It agrees reasonably 
well with the 3D secular theory for an inviscid disc. The 
slight offset of the precession frequency is attributable partly 
to errors in the secular approximation, which is valid only 
to leading order in e, and partly to the effects of viscos- 
ity. As described in Appendix the viscous damping of 
vertical motions considered in the full model can be rep- 
resented within the 3D secular theory by multiplying the 
coefficient 9/4 in the three-dimensional expression ^ for / 
by (l-i/3)/(l + i/3), where f3 = Qb + |a. Table[T]shows that 
this viscous secular theory gives good agreement with the 
full model for a = 0.1. 

The viscous damping rate of the modes is considerable. 
Although the dominant motion is horizontal and indepen- 
dent of z, and therefore does not incur any viscous forces 
in our approximation, the accompanying vertical motion is 
damped. To excite eccentric modes in a three-dimensional 
disc, this damping must be overcome. Viscous overstability 
may be able to do this, but detailed calculations are required 
and there is uncertainty in the applicability of a Navier- 
Stokes viscosity to turbulent stresses in the disc. A simple es- 
timate can be made as follows. In a two-dimensional, isother- 
mal, Keplerian shearing sheet with constant kinematic shear 
viscosity v and no bulk viscosity, the maximum local growth 
rate of the overstability is ~ 0.034 00 and occu rs for a ra- 
dial wavelength « 13ff (|Latter fc Ogilviell2006h . Ahhough 
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Table 1 . Scaled frequency eigenvalues obtained from the full lin- 
earized equations for a disc with e = 0.02, a = 0.1, Q = 20 and 
a = 2. Rapid convergence is seen with increasing values of the 
vertical truncation number A'^. Comparable results from the secu- 
lar theories are shown below. A negative value of lm{ui) represents 
the (scaled) exponential damping rate of the mode. 



Version 


Re((i) 


Im(a)) 


Full, N = 


1.189356 




Full, N = 2 


2.782638 


-0.449908 


Full, N = 4 


2.782585 


-0.449510 


Full, N = 6 


2.782585 


-0.449510 


2DS secular, inviscid 


0.636435 




2DO secular, inviscid 


1.193322 




3D secular, inviscid 


2.890299 




3D secular, viscous 


2.829523 


-0.448547 



our eigenfunctions do not have an obviously wavelike form, 
this suggests that overstability may be able to compensate 
for the damping rate found in Table [1] which corresponds 
to only ~ 0.0018 aQ for the parameters adopted there. 

A sufficiently large viscosity is required to couple differ- 
ent layers in the disc effectively. If the viscosity is reduced 
to Q = 0.01, with other parameters as in Table [T] a fre- 
quency eigenvalue of a; = 2.7831 — 0.1793 i is obtained, and a 
slightly larger value of N is required to obtain the same con- 
vergence. While the precession rate agrees well with the case 
of Q = 0.1, the damping rate is now larger than predicted 
by the viscous secular 3D theory. This happens because of 
the increasing z-dependence of the horizontal motion in the 
absence of a strong coupling between layers; although U2 is 
still much smaller than uo for a — 0.01, the ratio U2/U0 
is several times larger than in the case a — 0.1. The vis- 
cous damping of this shearing motion is more important, 
for a = 0.01, than that of the vertical motion considered in 
the viscous secular theory. 

If Q is reduced further, the precession rate starts 
to deviate from the 3D secular theory and more verti- 
cal structure develops in the eigenfunction, requiring a 
larger v alue of for convergen ce. Similar behaviour was 
found bv lLatter fc Ogilvid l|2006l ). The behaviour of a three- 
dimensional inviscid eccentric disc could be very difficult to 
describe. 

We also note that the stresses associated with tangled 
magnetic fields in a disc in which the magnetorotational in- 
stability occurs may provide an elastic, or v iscoelastic, cou- 
pling between different layers (|Ogilviell200ll ) . The associated 
damping rate may be smaller than that estimated on the ba- 
sis of a Navier-Stokes viscosity. 



5 CONCLUSIONS 

In this paper we have examined the linear dynamics of 
one-armed oscillation modes in the circumstellar discs of 
Be star s . A three-dimensional effect, first identified by 
lOgilvid (|200lD but neglected in previous treatments of 
Be stars, makes a crucial positive contribution to the preces- 
sion rates of such modes. It allows confined prograde modes 
to be obtained for all reasonable disc temperatures and stel- 
lar rotation rates. This property allows us to give a unified 



description of eccentric discs of Be stars of all stellar types 
without introducing uncertain radiative forces or modifying 
the inner boundary condition. While we do not deny that in 
some cases radiative forces might be important, or that the 
inner boundary condition might differ from a rigid one, we 
have shown that these innovations are unnecessary. We ob- 
tained these results using a secular theory of eccentric discs 
and confirmed them using a spectral treatment of the full 
linearized equations for three-dimensional isothermal discs 
including viscous terms that couple the horizontal motions 
at different altitudes. In order to make these modes grow, 
viscous damping must be overcome by an excitation mecha- 
nism such as viscous overstability, which will be investigated 
in a subsequent paper. 

The three-dimensional dynamics that we have described 
may also have important consequences for the behaviour 
of eccentric discs in other circumstances. For example, in 
cataclysmic variable stars exhibiting superhumps, the re- 
lation between the precession rate of the disc and the bi- 
nary mass ratio is an important observational property 
that is not fully explained by current theoretical models 
(Goodchild & Ogilvic 2006; Smith ct al. 2007). 

The physical model adopted in this paper is idealized. 
To explain in detail the observed cyclical behaviour of the 
emission lines of Be stars within this theoretical framework 
is likely to require a treatment of non-isothermal and non- 
linear effects as well as a better understanding of the time- 
dependent behaviour of the density distribution of the cir- 
cumstellar disc. In addition, future work should attempt to 
identify and assess mechanisms, such as viscous overstabil- 
ity, by which one-armed oscillation modes may be excited in 
circumstellar discs. Nevertheless, the effect investigated in 
this paper enhances the credibility of the one-armed oscil- 
lation model by providing a natural explanation of confined 
prograde modes. 
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APPENDIX A: DERIVATION OF THE 
SECULAR THEORY 

We consider equations (|48|l - (|51|l together with the viscous 
effects described in Section [4.21 When v„ and ui„ are elimi- 
nated, we have 



/ - 2 2\ 



2mn 



h„ 



liOhn rutin 

+ {dr\nH)[nhn + {n ^ l)(n + 2)hn+2\, 



hn rriK 

+ ^ „ . — u„ 



(Al) 



2rf2a)hr7 



+ — dr{rT.Un) + {drlnH){nu„ + u„^2) =0, (A2) 

where uihn ~ u + ianQ,z and livn = uj + i(ab + |a)(n — l)fiz 
are the Doppler-shifted wave frequency modified to allow 
for the viscous damping of horizontal and vertical motions, 
respectively. We are interested in the case m — 1. 
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This system of equations is not closed because the uq 
equation refers to /i2, which depends on U2. In turn, the U2 
equation refers to /14, which depends on U4, and so on. How- 
ever, suppose for the time being that U2 can be neglected 
in the /12 equation (i.e. equation IA2I with n — 2). Then we 
have a closed system 



— -uo = —drho H -ho + 2(9r ln_H")/i2, 



4i) 


ho 




ci +2 




cl 1 


iiiv2 


iiih2 



tUo + —^dr{rT.uo) = 0, 



^ + {dr\nH)uo 



0. 



(A3) 
(A4) 

(A5) 



These equations agree with a two-dimensional theory except 
for the additional /12 term in equation (|A3|I . which can be 
related to uo through the approximate equation (|A5[) . 

In the secular approximation we neglect the differences 
between f2, k, and Q.k except where essential, so that 
Q. oc r-^/'^ and H oc r ' , and we assume a low frequency 
\uli\ <^ SI. The leading approximations to the above equations 
in the inviscid case are then 



2i 



2n 



ui uo 



-drho 



2ho 



3h2 

r 



= 0. 



(A6) 



(A7) 



(A8) 



These combine to give equation © for u — uq, with the 
three-dimensional expression ([8]) for /. If h2 is neglected 
altogether, equation © is obtained, but with the two- 
dimensional expression ([7} for /. It can be seen from the 
above equations that /12 always makes a positive contribu- 
tion to the precession rate uu. 

Neglecting U2 compared to uo is equivalent to assuming 
that the eccentricity is independent of z, since Heo(C) ~ 1 
and IIe2(C) = C'^ ~ 1- Under what conditions is this assump- 
tion reasonable? Rough estimates based on equations (|Aip 
and (|A2|) show that 112 can be neglected in the /i2 equa- 
tion if ail is much larger than the precession frequencies 
u> or /, meaning that the shear viscosity prevents the sig- 
nificant development of a z-dependent horizontal motion. 
This condition is readily satisfied in the circumstellar discs 
of Be stars for reasonable values of a. Ultimately, however, 
the validation of this approximation comes from the numer- 
ical solution of the full system of linearized equations. 

If viscosity is retained in this analysis, within the secular 
approximation, then the coefficient of /12 in equation (lASP is 
multiplied by (1 -I- i/3)/(l — i/3), where 13 = Ob + |a. In this 
case, equation (|6} is again obtained, but in expression ([8]) 
for / the coefficient 9/4 of the three-dimensional term is 
multiplied by (1 — i/3) /{! + i/3). In this case the mode decays 
because of the viscous damping of the vertical motion. For 
sufficiently small a, however, the decay rate is enhanced by 
the viscous damping of the z-dependent horizontal motion. 
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